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Abstract 

Low accuracy of the Solvation Free Energy (SEE) calculation is a known problem of the 
numerical methods of the Integral Equation Theory of Liquids and the Classical Density 
Eunctional Theory (Classical DET). Although functionals with empirical corrections can es¬ 
sentially improve the predictability of the methods, their universality is still a question. In 
our recent paper we connected the SEE calculation errors with the incorrect pressure in the 
Classical DET and proposed the a posteriory correction to improve the results (J. Phys. 
Chem. Lett., 5, 1925-1942 ). This paper raised a discussion in the community. In particular, 
recently appeared a critical reply where pointed some thermodynamical inconsistencies of the 
derivations in our paper (J. Chem. Theory Comput., 11, 378-380). In the present work we 
re-derive the pressure correction in a more simple way and show that despite the inaccuracies 
during the derivation, the final form of the previously derived correction is correct. We also 
test the applicability of the proposed correction to the functionals which include a three- and 
many- body terms from the fundamental measure theory (EMT) for hard sphere fluid. We 
test all the functionals on a set of model systems and discuss the obtained results. 


1 Introduction 

Solvation Free energy fundamental quantity in physical chemistry which allows ones to compute 
many useful properties of substances in solution Unfortunately, the solvation free energy calcula¬ 
tions are computationally expensive. That’s why there exist alternative, faster but sometimes less 
accurate methods, like continuum electrostatics models[T] or energy representation technique[2]. 
In this group of methods, the methods based on the Classical Density Functional Theory (DFT) 
are one of the most promising techniques for solvation free energy calculations. They are orders of 
magnitude faster than the simulations bu nevertheless preserve an important information about 
the solvent structure around the solute [3]. 

However, the accurate solvation free energy calculation in the Classical DFT approach appeared 
to be a challenging task. Although the methods based on the simplest Hyper Netted Chain (HNC) 
approximation lai can give a qualitatively correct estimation of the solvent structure around the 
solute the calculated values of the solvation free energies are unrealistic and far from experimental 
ones |S]. It was noticed, that the error in the solvation free energy calculations is proportional to 
the partial molar volume of the solute[6]. There were proposed the empirical correction models for 
the RISM and 3DRISM, which used this observation and were able to predict the solvation free 
energies of simple solute with accuracy of 1 kcal/mol[71 El |9]. In our recent paper we were able 
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to show that these empirical corrections can be explained by the difference in pressure in HNC 
approximation and in experiments [10] . 

Another approach to the solvation free energy calculations is to improve the HNC model itself. 
The aim of improvements is to include the three and many body interactions to the model. Some 
of the models include the empirical terms for improvement the water structure [iiiiniDS], another 
use the three and many-body interactions from the simplified models . The models which use the 
experimental c-functions, and include the many-body interactions from the Fundamental Messure 
Theory (FMT) are called FMSA or the models with the Hard Sphere bridge [TTl 115] . It was 
shown, that using these models it is possible to predict accurately the solvation energies of alknes 
[IS], and also more complex solutes im. One of the disadvantages of these models is the fact 
that they still contain an empirical parameter: hard sphere diameter of the solvent. The result of 
the calculations are very sensitive to this parameter and the optimal choice of the proper solvent 
diameter is a challenging task. Another disadvantage of this kind of empirical methods is that the 
predictability of the model with fitting parameters is always a question, and need to be proved 
by the tests on extensive number of systems, which is computationally demanding and never give 
a full understanding of the model limitations. In a light of our previous paper it would be also 
interesting to measure the pressure in these models, to answer the question whether it is realistic 
or does the HS diameter is just another fitting parameter. 

Although the explanation given in our previous paper is not perfect and we agree with the 
derivations given in [18] we still cannot agree with the conclusion of that paper that the proposed 
correction does not have a physical grounds. We show below that the final formula for the solvation 
free energy calculation given in [10] can be correctly explained by a simple difference in pressure 
in the HNC model (which is extremely high) and atmospheric pressure which is almost negligible 
in comparison to it. So, the numerical results of the paper [10] are still correct and transferable to 
any Classical DFT calculations with HNC functional. 

Additionally, we derive the similar pressure expression for the functionals with the HS bridge 
and demonstrate that the sensitivity of the model to the hard sphere radius is also due to the 
pressure of the model. We discuss the different ways to choose the solvent hard sphere diameter 
which can either reproduce the correct pressure or the correct free energies of the small solutes 
(with the price of loosing the physical consistency). At the very end we try to find the systematic 
components in the expressions which reproduce the correct pressure, and show that for huge variety 
of systems it can be explained with a good accuracy by adding the ideal-pressure correction, which 
poses new challenges for the theoretical justification of this correction. 

2 Pressure correction in classical DFT and related meth¬ 
ods 

Let us consider the process of solvation of the rigid solute. As it was absolutely correctly mentioned 
in Ref. [18], it does not matter which kind of process (NVT, NPT, fiVT) do we consider, we will 
always recover with the same solvation free energy. So, let us consider the solvation process in 
NPT ensemble, where the pressure, temperature and the number of particles are constant as this 
ensemble is commonly used in experiments and simulations. By definition in the solvation process 
the solute is transfered from the fixed position in the ideal gas to the fixed position in the solvent 
[T9] . Then during the solvation process the volume of the liquid system changes by some volume 
AH, which is (by definition) the partial molar volume of the solvent (see Figure [^. According to 
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Figure 1: Solvation process in the NPT ensemble and three contributions to the Gibbs energy 
change. 
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Figure 2: Radial distribution functions of water around methane computed with different Hard 
Sphere radius parameter. 
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the basic thermodynamic relations we can hnd the components of the Gibbs solvation energy AG: 

AG = AF + PAV (1) 

where AF = AU—TAS is the Helmholtz free energy change, P is the pressure. Different properties 
of the liquid contribute to the AF and PAV components of the solvation free energy. Both: 
internal energy change AU and the entropy change AS depend in the first turn on the structural 
changes in the solvent, but not directly depend on the volume change AV. In turn, the pressure 
term PAV is responsible for the energy change due to the volume change, while it does not depend 
on the structure of the solvent. Separation of the energy contributions in this way allows us to 
propose the reasonable approximations for the solvation free energy AG. 

It is known, that the methods based on the Classical DFT approach can qualitatively correctly 
predict the solvent structure around the solute [2D1 ED ED ED Elj- For example, in Figure 
are shown the radial distribution functions of water around the methane in water. Although 
the quantitative behavior is not so good, and there is a series of well-known problems of the 
Classical DFT methods, we may still assume that the the structure predicted by the Classical 
DFT methods is good enough to roughly reproduce the structure-based component AF of the 
Solvation free energy. 

In the other hand, it is also known and shown in our previous paper [10] , that the models based 
on the Homogeneous Reference Fluid approximation, like HNC models, fail to predict the correct 
pressure. As we show below, the typical pressure in HNC approximation with SPCE solvent at 
room temperature and usual water density is about 11.5 KBar, which is 4 order of magnitude higher 
than in experiments. In this situation it is reasonable to expect that the pressure contribution is a 
dominant contribution to the solvation free energy and the main source of the errors of the model. 

To correct these errors we can have several approaches. The simplest one is an a posteriori 
pressure correction. Writing equation Q for both: DFT calculations and experiments and using 
the approximation AF^pj’ ^ AF^^pt we get the approximation for the experimental solvation free 
energy AGexpt- 

AGexpt ~ AG DFT — PdftAVdft + PexptAV^xpt ( 2 ) 

Another approach is to use the advanced free-energy functional. As the functionals are rarely 
completely parameter-free, we have to possibilities: either to £t the parameters to get the correct 
SFE of some given molecule (set of molecules), or take the parameters which reproduce the correct 
pressure. (The ideal functional would do both.) 

In the next sections we consider the described approaches and discuss the numerical results. 


2.1 HNC and a posteriori Pressure Corrected HNC 


One of the simplest and most popular models considered in the classical DFT theories is the Hypper 
Netted Chain (HNC) or Homogeneous Reference Fluid (HRF) approximation. In this model the 
excess free energy functional is expressed as a second-order Taylor series approximation around 
the homogeneous fluid density po- Subtraction of the piN term with p = and reference 

fluid free energy leads to the following Grand potential functional with respect to the reference 
fluid: 


AD[p] = D[p]—D[po] = kT 


p(l)lndh_Ap(l) 

Po 


/ kT 


Ap(l)c(12)Ap(2)dld2 


where arguments 1, 2 stay for the positions and orientations of the solvent molecules, U is the 
potential of the solute molecule, Ap(l) = p(l) — po, c(12) = — (36"^ / 6p{l)6p{2) is a pair direct 
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Figure 3: Results of the MDFT calculations for the small Alkanes and rare gases. HNC corre¬ 
sponds to the Homogeneous Reference Fluid approximation, HNC/PC is the same but includes the 
pressure correction, HS/FlT corresponds to the calculations with the Hard Sphere bridge, where 
the Hard Sphere radius is chosen to reproduce the SFE of methane, HS/PRES is the same, but the 
Hard Sphere radius is chosen to reproduce the Atmospheric pressure, HNC/PC-I- and HS/PRES-I- 
are the results with additional pkTAV term. 

correlation function, A; is a Boltzmann constant, T is a temperature, (3 = 1/kT. By minimizing 
the functional over the solvent density p(l) one can find both: the solvation free energy AH and 
the density distribution p(l). 

It is known that without the pressure correction HNC calculations dramatically overestimate 
the solvation free energies of the solutes EE]. To demostrate that fact we performed the HNC 
calculations for small hydrophobic solutes in water (see Figuresquares). As expected, the HNC 
SFE predictions are always too high, sometimes by 200-400%. 

To check whether the errors are random or systematic, we performed the calculations for the 
simpler system: hard sphere solutes in water. It is known, that for the macroscopic solutes the 
free energy is proportional to the surface of the solute. In case of the liquid-vapour interface the 
proportionality coefficient is the surface tension 7 . In Figure we see the dependency of the 
Solvation Free Energy per surface area on the hard sphere diameter. We see that HNC functional 
fails to produce the correct behavior of the surface tension: the solvation free energy per surface 
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area in HNC approximation is proportional to the hard sphere radius and never comes to the 
plateau. So, HNC fails to reproduce even the qualitative behavior of the solvation free energy for 
the simple model system. 

To reduce the errors, we can use the a posteriori pressure correction described in the previous 
section. The pressure in the HNC approximation can be found using the known thermodynamic 
relation for the homogeneous fluid: 

H[po] = -PV 

If we set zero density to the functional, we get: 

kT 

Phnc = AH[0]/H = -n[po]/V = pokT - —plc{k = 0) (3) 


where c{k = 0) = f c(12)dl. The definition of the pressure, given above corresponds to the 
results in Refs. [25l 1^ . We should say, that this is notthe only possible definition. To avoid 
misunderstanding in the community, we briefly discuss the different ways of pressure dehnition in 
Classical DFT and their domain of applicability. 

The pressure computed in ([^ is a so-called compressibility route pressure, which differs from the 
virial route pressure. We emphasize, that it is consistent in our derivations to use compressibility, 
but not virial route pressure. To support that, we can write the virial-route pressure in the simplest 
case of spherical particles [27] : 


P = pokT - 



du{r) 

dr 


g{r)dr 


where u{r) is a pair-wise potential and g{r) is a radial distribution function of the solvent. It is 
well seen that this quantity depends on the structural properties of the system only, and is the 
same for any choice of the approximate functional Ar2[p]. Thus, it cannot be used to correct 
its deficiencies. Also, the pressure in (|^ should not be mixed with the pressure which can be 
computed from solving the HNC equations for solvent, as it is done in Ref. [28]. In our model the 
solvent direct correlation function c(r) is extracted from the simulations, so it is inconsistent to 
use it in any way in the HNC equations for solvent. In turn, the pressure defined as ([^ correctly 
defines the energetic behavior of the system, in a sense that it shows the energy change then the 
volume of the system changes, and this is exactly the quantity we are interesting in. 

For the SPCE water model used in our simulations the compressibility route pressure computed 
by ([^ we get the the pressure of 11.477 KBar. As this pressure is 4 orders of magnitude higher than 
ones used in experiments, we can neglect the PgxptAWxpt term in (|^. In that case the solvation 
free energy in the HNC approximation with the pressure correction (HNC/PC) can be written as 


^Ghnc/pc — Ar2[p] — pQkT{l — ^c{k — 0 )AVhnc 


This result is completely consistent with the result which we got in Ref [10], (eq (11)), which 
unfortunately had not the best reasoning. 

Looking at Figure]^ we see, that the results for HNC/PC improves the solvation free energies 
of small solutes (rare gases and methane), but it essentially underestimates the solvation energies 
of larger Alkanes. Also, as it can be seen in Figure]^ the HNC/PC model is able to reproduce the 
qualitative behavior of the free energy of a big hard sphere solute in water. 

Despite the fact that pressure correction can to some extent correct the misleadingly high 
error of the HNC model, it is still far from being perfect. The assumption used in the pressure 
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Figure 4: Dependence of the free energy per surface area for the Hard Spheres of different diameter 
immerced in water. Results are shown for different free energy functionals. The abbreviations are 
the same as in Figure 


correction approximation is that HNC can give a relatively good approximation of the solvent 
structure around the solute. Although in many cases this assumption is quite reasonable, in some 
situations the structure in HNC approximation is very far from reality. One of such cases is the 
contact density of the water with a hard-sphere. According to the so-called contact value theorem 
the contact density near the hard wall is proportional to the pressure [T]. Note, that for the 
compressibility route pressure and approximate functional it should not necessarily hold, but, a 
high contact density can be another indication of a functional failure to reproduce the simulation 
and experimental results. 

In Figure the contact density for the different hard sphere diameters is presented. We see 
that the HNC density grows up to g{rcontact) = 7-7- The real contact density should decrease 
with the hard sphere diameter increase, and tend to the small value which corresponds to the 
atmospheric pressure P/kT cs 0.01. The water RDF for the Rbaii = 12 A is shown also in Figure 
We see that the radial distribution function (RDF) calculated with the HNC functional shows 
completely wrong behavior near the hard sphere surface. So, we may expect also a big errors in 
calculation the thermodynamic quantities as well. So, to conclude: it is known and proven many 
times that HNC functional is not perfect. It predicts the solvation free energies proportional to 


7 










10 


u 

-p 

c 

o 

u 

ip 

Cn 


HNC 

HS/FIT 

HS/PRES 


8 - 


6 - 


4 - 


2 - 



-e o o e 


10 


20 30 


40 


50 


60 


70 


80 


^ball 


Figure 5: Contact density and the maximum density as a function of hard sphere radius for 
different free energy functionals. Note: the contact density coincides with the maximum density 
for HNC and HS/FIT functionals. The abbreviations are the same as in Figure]^ 


the volume, and thus fails to reproduce the surface tension. Although the pressure corection can 
correct the unnaturally high values and make them more reasonable, still the deficiencies of the 
HNC functional do not allow to improve the method further. HNC fails to reproduce the correct 
RDF, and in many cases the difference is drastic, so we cannot expect the correct thermodynamics. 

2.2 Hard Sphere Bridge functional with empirical solvent diameter 

As it was discussed above, there is an alternative way to improve the solvation free energy calcula¬ 
tions. To do it one can use an improved functional, which includes many-body terms. In our paper 
we use the mixed functional approach: use the two-body c-function extracted from the simulations 
but add the many-body terms from the model where the analytical solution is known. The good 
candidate is the fundamental measure theory (FMT) with the functional for the hard sphere fluid. 
In that case we dehne the full free energy functional in the following way: 
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Figure 6: Density distribution functions of water around the Hard Sphere of radius R=12 A , 
calculated with three different functionals. The functional name abbreviations are the same as in 
Figure 


where ^ ^ many body interactions of hard sphere fluid, which includes terms starting from 

the 3’’*^ order and higher: 


r -pexc r r r 

J^Hp[p\= J^Hs[p]-J^Hs[Po]-—^ J Ap{r)dr+— JJ Ap{ri)cHs{\r 2 -ri\)Ap{r 2 )dridr 2 (5) 


CHs{\r2 - ril) = 


£2 -pexc 


Sp(ri)Sp(r2) 


( 6 ) 


PO 


In the paper we use the Percus Yevick functional [2] in the scalar representation 


-pexc _ 

^HS — 



(7) 
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Figure 7: Free energy of methane calculated with Hard Sphere bridge functional with different 
hard sphere radii (HS). HS + PC - the same values with the Pressure Correction, where the 
pressure contribution to the free energy is subtracted. 

where Uj(r) = fy p(r')cji(r — r')dr' and spherically symmetric weight functions w(r) are dehned in 
the fc-space as follows: 


kR 

oJo{k) = cos kR + sin kR 

(9) 

ui{k) = —{sin kR + cos kR) 

(10) 

^ 2 {k) = — smfcit 

k 

(11) 

Att 

Cj 2 ,{k) = — {sin kR — kR cos kR) 
k'^ 

(12) 


where R is the radius of the hard sphere. 

The hard sphere radius is not pre-dehned for the aqueous solution, so we can consider it to be 
a free parameter. In Figure the dependence of the Solvation Free energy of methane on R is 
shown. We see that free energy of solvation depends on the hard-sphere radius R, especially in the 
region 1.2A < R < l.sA where solvation free energy is extremely sensitive to it and change by 40 
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kJ/mol. So, we can choose R which reproduces the solvation free energy of methane. In our case 
it is i? ~ 1.41A. If the hx this parameter, we obtain the functional which can be used to calculate 
solvation energies of arbitrary molecules. We refer this functional as HS/FIT. 

Looking at Figure we see, that HS/FIT functional is able to reproduce the MD solvation 
free energies with a good accuracy. It is almost perfect for the small one-atom molecules. For 
larger Alkanes there is a growing deviation, which is still much smaller than HNC/PC results. In 
principle, as it was demonstrated by Wu 021 , using this method with a proper parameterization 
it is possible to reproduce the SFE for a larger sets of molecules with a quite high accuracy. In 
Figure|^we see that HS/FIT method reproduces the height of the hrst peak of methane-water RDF 
almost correctly. However, it is too shallow, and the position of the hrst minimum is incorrect. 
But in general we may conclude that the results for the small solutes are much better than HNC 
and nearly of the accuracy of simulations. 

However, result for the Large Hard Sphere Solute in water are different. In Figure]^ we see that 
the HS/FIT functional does not reproduce the correct solvent structure around the hard sphere 
solute. In Figure we see, that the HS/FIT functional gives the contact density smaller than 
HNC, but still too high to correspond to real pressure. The fact that the HS/FIT pressure is 
incorrect is even better seen in Figure where we see that HS/FIT the solvation free energy is 
proportional to the volume of the solute. 

We can conclude that using the the hard-sphere bridge functional with a radius parameterized 
over the SFE of small molecule (in our case methane), we can recover the SFEs of other small 
molecules (noble gases and Alkanes) with the good accuracy. However, in this way we are unable 
to recover the correct pressure, which results in the incorrect Solvation Free energy per surface 
behavior and incorrect solvent structure for the big hard sphere solutes. So, also HS/FIT method 
can be used to reproduce the SFEs of small solutes, one need to remember about its limitations 
of applicability. 

2.3 Hard Sphere bridge functional with correct pressure 

Instead of using the HS/FIT functional, we can also construct the functional which reproduces the 
correct pressure of the system. To do it, we hnd the compressibility-way pressure produced by the 
functional (|^. As in case of HNC functional, we use the relation Af 2 ji/ 5 [ 0 ] = PH, which gives the 
following pressur^ : 

kT 

P = pokT - —pl{c{k = 0) - CHsik = 0)) + P“^^ 

where 

PpMT = PokT{AB + QB^ + 3P3) 
where B = pW^/{l — pHs), H 3 = 4 / 37 rP^. 

This allows us to plot the pressure in HS-bridge functional as a function of R. In Figure we 
see that dependency. We see that for R < lA the pressure is almost the same as in HNC model, 
however in region I. 2 A < R < l.sA the dependency is very strong. We see that the pressure 
in HS/FIT model with R = 1.41 is 4.56 KBar. The radius R = 1.4574A corresponds to the 
atmospheric pressure. We refer the hard sphere bridge functional with with radius as HS/PRES. 

To check whether our derivations are correct and the HS/PRES functional really reproduces 
the correct pressure of the system, we can look at the Figure We see that at large Hard 

^See Supporting Information for the derivation 
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Figure 8: Pressure for the functional with hard sphere bridge as a function of the hard sphere 
radius of FMT functional. The pressure of HS/FIT functional with i? = 1.41 A is 4.56 KBar. The 
pressure of 1 Bar corresponds to Rhs = 1-4574 A . 

Sphere diameter the HS/PRES free energy is proportional to the surface area, but not to the 
volume, which is an expected result. Now we can check how good can the HS/PRES functional 
reproduce the solvent structure. In Figure]^ we see that the hrst peak of methane-water RDF is 
not reproduced completely correct: it is shifted from the methane and is a little bit smaller than 
in simulations. However, the hrst minimum and the second peak are reproduced much better than 
in other methods (HNC, HS/FIT). The results for the large hard spheres in water are even better. 
In Figure]^ we see the correct behavior for the contact density for HS/PRES functional: it tends 
to a small value when the radius of the hard sphere increases. In Figure|^we see, that Hard-sphere 
water RDF for HS/PRES methods behaves qualitatively correct and much more consistent with 
the simulations than other methods. However, solvation free energy results are not so good. In 
Figure we see that the HS/PRES functional underestimates the free energies of all solutes, it is 
especially visible for multi-atom Alkanes. 

As the pressure correction ([^ can be dehned for any functional, we in principle, we can also 
imagine a hybrid method, which uses both: the pressure correction and the hard sphere bridge. In 
Figure we see that this functional underestimates the SFE for any reasonable parameter R. The 
only way to recover the correct solvation free energy is to use the R > lA, which corresponds to a 
negative pressure and is thus completely unphysical. Unfortunately, with the hard sphere bridge 
functional we cannot get all: good solvation free energy, correct structure and thermodynamically 
consistent pressure. However, we may decide what is more important for us in each concrete case 
and use the corresponding parameters. 
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Figure 9: Dependencies of the errors of the functionals after the pressure correction on the partial 
molar volume of the molecules. HNC/PC+ and HS/PRES+ corresponds to the calculations with 
the additional pkTAV term. 

3 Additional correction 

In the previous section after analysing the ability of different functionlas to reproduce SFE and 
testing them on the model system of large hard spheres solvated in water we came to the conclusion 
that we cannot simultaneously accruratly reproduce all the properties: pressure, SFE and solvent 
structure, and thus need to look for a compromise. However, we need to say that this conclusion is 
not completely right. Doing the numerical experiments for the systems with different parameters 
of the functionals we noted that the errors in SFE calculation of the functionals which correctly 
reproduce the pressure (HNC/PC, HS/PRES) is not random. 

In Figure one sees the dependency of this error on the partial molar volume AV. We can 
spot an evident linear dependency for both of the methods. The slope of the dependency is almost 
equal for both of the methods and can be approximated with a good accuracy by the expression 
pkT. Doing the same analysis for the greater number of molecules we can suppose that this is not 
just a coincidence, but a systematic behavior (See Supporting Information). On the basis of our 
considerations we can propose an additional correction: 

AGpc+ = AD[p] - (P - Pe.pt)AV + pkTAV (13) 

This correction can be used for any functional. We denote the results where this correction was 
used by writing an additional “+” at the end, for example HNC/PC+, HS/PRES+. In Figure 
we see that the additionally corrected results are much less dependent on AV and are much more 
close to the experimental ones. 
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To check the effectiveness of the additional correction we performed the calculations for the 
larger set of molecules. In Figure Figure ^we see comparison of the different methods with 
Pressure Correction for 122 organic molecules ^ . We see that results with the additional correction 
can essentially improve results for both cases: a posteriori pressure corrected HNC functional 
(HNC/PC+) and HS-Bridge functional with the correct pressure (HS/PRES+). In the case of 
HNC/PC+ the resulting error is comparable to the error of simulations. We see that in general 
HS/PRES+ gives a more dispersed results than HNC/PC+. However, this can be explained by 
the numerical instability of the HS-bridge calculations for strongly polar molecules. For non-polar 
molecules with possitive SFE results of the both functionals are comparable and have a similar 
accuracy. So, the results suggest us that the correction (13) is not artificial but possibly have some 
well physical grounds and can be successfully used in Classical DFT calculations. 

We need to comment separately the recent paper Ref [30], where the Pressure Correction 
method is used with the 3DRISM calculations. In that paper two variants of the a posteriori 
pressure correction: with and without additional pkT term ( which corresponds to HNC/PC and 
HNC/PC+ in our definitions). The numerical results presented in this paper suggests that the 
HNC/PC produces more reliable results in the 3DRISM case. As an explanation of this fact it 
was stated, that 3DRISM is forulated for the Grand Canonical ensemble (in contrast to MDFT, 
as it could be understood from the context). Although we do not deny the computational results 
of Misin & coathours, we need to note, that the explanation of the difference in MDFT and 
3DRISM results, given in the above mentioned paper is not good enough. In particular, for the 
case of spherical solvent 3DRISM equations are completely equivalent to the MDFT and can be 
obtained from the later by taking the functional derivative and equating to zero the derivatives 
over density function. So, at least for this case no differences in results are possible. From our 
point of view it is much more reasonable to assume that in the 3DRISM model the pressure differs 
from the MDFT pressure. This assumption seems even more reasonable, if we take into account 
that 3DRISM theory describes the mixture of sites connected in the complicated way, so it is 
completely possible that each site can give its own contribution to the pressure. In such a way, we 
argue that the conclusions of Ref [Misin] do not contradict our results for the MDFT functional 
and only additionally indicate that the question of the applicability of the Pressure Correction to 
the 3DRISM calculations should be the topic of additional investigation. 


4 Conclusion 

In our paper we performed an efficiency test of different Classical DFT functionals for the Solvation 
Free Energy calculations. We confirmed the known that that the HNC approximation gives a 
hugely overestimated predictions of the SFE of the small molecules. Using the model system of 
hard sphere solvated in the water we demonstrated that HNC functional cannot even qualitatively 
correct predict the dependency of the SFE on the size of the Hard Sphere. We have shown that 
for the large spheres the SFE calculated with the HNC is proportional to the volume, instead 
of the surface area of the sphere. We suggested that this can be connected with the fact that 
HNC approximation incorrectly reproduce the pressure of the system and as a consequence of it 
the —PAV contribution to the total solvation energy. We derived the expression for calculation a 
compressibility-route pressure in the HNC model, and have shown that for the SPCE water used in 
our calculations it is 11.5 KBar. Thus the PAV term is dominant in the error of SFE calculations. 
To improve predictions of the solvation free energies we proposed to use an aposteriori pressure 

^Experimental and simulation data taken from |29| . See the full list of molecules in Supporting Information. 
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Figure 10: Error of the free energy calculations for 122 molecules calculated with the free energy 
functional with the hard sphere bridge. Blue diamonds (HS/FIT) are the results for the Rhs = 
1.41 A which reproduces the free energy of methane. Green triangles (HS/PRES) are the results 
for the Rhs = 1-4574 A , which reproduce the correct pressure. Red squared (HS/PRES+) are 
the results which reproduce the correct pressure with additional pkTAV contribution added. 


correction (|^. We showed that HNC functional with the Pressure Correction (HNC/PC) correctly 
reproduces qualitative behavior of the SFE of hard spheres in water. However, this correction is 
unable to overcome all the defficiencies of the HNC model. In particular we demonstrated that 
the solvent structure near the hard sphere is not reproduced correctly which also should have 
effect on the SFE calculations. Although the errors of the SFE predictions of the small molecules 
in HNC/PC are much smaller than in HNC approximation, nevertheless they are still essentially 
large. To improve the HNC/PC functional we used the Hard Sphere Bridge (HS-Bridge) functional 
Q which includes the many-body correlation functions of hard spheres derived from the FMT 
theory. This functional depends on the parameter R - hard sphere radius of the solvent molecules. 
We parameterized the functional to reproduce correctly the SFE of methane in water. We showed 
that using this parameterization allows one to reproduce correctly SFEs of the small molecules in 
water. However, this functional incorrectly reproduces the solvent pressure and as a result gives 
an incorrect behavior for the SFEs of large hard spheres in water. This functional is also unable 
to reproduce the distribution of the solvent molecules near the hard sphere ball. 
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Figure 11: Error of the free energy calculations for 122 molecules calculated with the Homogeneous 
Reference Fluid (HNC) free energy functional. Orange circelss (HNC/PC) are the results for Ho¬ 
mogeneous Refenrence Fluid approximation with the Pressure Correction. Violet Stars (HNC/PC 
-|- rhokT) are the same result with additional pkTAV contribution. Green pluses (MD) are the 
results of the MD simulations. 

We used an alternative approach and performed another parameterization of the HS-Bridge 
functional, to get the correct pressure of the system. The functional parameterized in this way 
can much better reproduce solvation parameters of the hard sphere in water. Both: solvation 
free energy and solvent structure behavior are qualitatively correct. However in this case SFEs 
of the small molecules in water are essentially underestimated. We additionally tested “hybrid” 
functionals, which include both: HS-Bridge and Pressure Corrected and showed that for any 
reasonable value of the R-parameter these functionals always underestimate the solvation free 
energy of small solutes. In such a way we came to the conclusion that we need to look for a 
compromise between the thermodynamical consistency and relevant results for the model systems 
and accuracy of the SFE predictions for the “real” systems like hydrated organic molecules. In 
addition, we studied dependency of the errors of the SFE predictions with the functionals which 
reproduce the correct pressure on the partial molar volume of the molecule, and noted that the 
errors can be reasonably well corrected by adding additional pkTAV term. On the basis of these 
considerations we proposed an additional correction which includes that component. To check 
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whether this correction is universal or not we performed the SFE calculations for 122 organic 
molecules, and showed that for both studied functionals results with the additional correction are 
the best. In case of the HNC/PC+ functional the errors of predictions are comparable with the 
errors of simulations. So, this allows us to suggest that this correction should be used in the MDFT 
calculations to achieve the best accuracy of the SFE predictions. The choice of the functional - 
HNC/PC+, HS/PRES+ or other depends on the particular task, and is the matter of question 
which of the thermodynamical properties are important in the particular research. 
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